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PROBLEM  STATEMENT 


This  research  project  investigated  whether  hydrologic  model-based  interpretation  of 
transient  GPR  (ground-penetrating  radar)  data,  the  ground  wave  in  particular,  can  lead  to 
improved  characterization  of  soil  water  dynamics.  The  strong  dependence  of  the  high-frequency 
dielectric  permittivity  of  a  soil  on  water  content  has  made  electromagnetic  measurements  (e.g., 
time  domain  reflectometry  (TDR),  ground-penetrating  radar  (GPR),  and  satellite  radars)  popular 
tools  for  monitoring  near  surface  moisture  conditions  [e.g.,  Huisman  et  ah,  2001;  Entekhabi  et 
ah,  2004].  The  GPR  ground  wave  is  of  particular  interest  for  providing  high-resolution  maps  of 
near-surface  water  content  at  the  catchment  scale  [Grote  et  ah,  2003;  Weihermuller  et  ah,  2007]; 
such  data  would  be  invaluable  for  applications  including  understanding  local  rainfall-runoff  and 
infiltration  processes  as  well  as  calibrating  satellite-based  radar  measurements  for  real-time 
monitoring  of  regional  moisture  conditions.  A  limitation  of  current  interpretation  techniques, 
however,  is  the  assumption  that  the  effective  permittivity  derived  from  a  GPR  measurement  can 
be  used  to  estimate  the  average  water  content  of  a  soil,  which  can  subsequently  be  used  to 
interpret  hydrologic  processes.  Previous  studies  have  shown  that  problems  like  preferential 
sampling  of  high  velocity  zones  by  EM  waves  [Galagedara  et  ah,  2005a,b;  Moysey  and  Knight, 
2004]  and  dispersive  wave  guides  resulting  from  the  infiltration  of  water  [e.g.,  van  der  Kruk, 
2006]  can  lead  to  systematic  biases  in  water  content  estimates.  This  problem  is  compounded  if 
the  biased  water  contents  are  then  used  to  calibrate  non-linear  unsaturated  flow  models. 
Following  this  sequential  approach  to  geophysical  data  integration,  i.e.,  estimating  transient 
water  contents  with  GPR  and  subsequently  performing  hydrologic  analysis,  may  therefore  lead 
to  both  inaccurate  estimates  of  water  content  and  poor  predictions  of  unsaturated  flow  and 
transport.  This  research  used  modeling  studies  to  evaluate  whether  a  new,  coupled  approach  to 
hydrogeophysical  estimation  of  hydrologic  parameters  could  outperform  the  sequential  approach. 
The  modeling  was  complimented  by  laboratory  and  field  studies  to  assess  whether  the  coupled 
approach  could  also  be  applied  to  real  data. 


SIGNIFICANCE  OF  THIS  RESEARCH  TO  THE  ARMY 

Developing  a  methodology  for  the  hydrologic  analysis  of  GPR  data  to  characterize  soil  moisture 
dynamics  is  of  critical  importance  to  the  Army  for  several  reasons,  including  understanding  risks 
related  to  subsurface  transport  of  Army-related  contaminants,  evaluating  ground  conditions 
affecting  troop  maneuvers  under  changing  environmental  conditions,  and  predicting  surface 
runoff  generation  during  rain  events  that  can  lead  to  flooding.  Understanding  the  response  of 
ground  penetrating  radar  under  dynamic  soil  moisture  conditions  is  also  of  the  utmost 
importance  for  improving  the  Army’s  ability  to  detect  targets,  such  as  landmines  and  buried 
ordinance,  under  field  operating  conditions. 
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SUMMARY  OF  MAJOR  RESULTS 

This  summary  provides  a  brief  overview  of  work  completed  and  results  achieved  in  this 
research.  A  more  detailed  review  of  methods  and  analysis  are  given  in  the  technical  report  in  the 
appendix. 

A.  Numerical  Comparison  of  Sequential  and  Coupled  Inversion  Schemes 

The  first  problem  addressed  in  this  research  was  to  assess  the  ability  of  sequential  versus 
coupled  inversion  schemes  to  constrain  the  parameters  of  an  infiltration  model  using  GPR 
traveltime  data.  Analytical  models  were  selected  for  this  analysis  to  allow  for  a  Markov  Chain 
Monte  Carlo  (MCMC)  analysis  to  be  performed.  The  Philip  [1957]  infiltration  model  was 
selected  with  a  rectangular  drainage  model  for  redistribution,  which  results  in  5  unknown 
hydrologic  parameters  controlling  flow  (sorptivity  (S),  saturated  hydraulic  conductivity  (Kg), 
initial  water  content  (0i),  saturated  water  content  (0s),  shape  parameter  for  K(0)  function  (N)). 
This  hydrologic  model  produces  a  sharp  wetting  front,  which  allows  for  calculation  of  GPR 
traveltimes  (i.e.,  air  wave,  ground  wave,  reflected  wave,  refracted  wave,  and  multiples)  directly 
using  ray  theory.  The  MCMC  analysis  was  performed  using  the  Metropolis-Hastings  algorithm 
to  sample  the  posterior  distribution  of  model  parameters  subject  to  the  GPR  data.  In  the  case  of 
sequential  inversion  GPR  traveltime  data  were  used  to  estimate  the  depth  of  the  wetting  front  and 
the  water  content  above  and  below  the  wetting  front,  which  were  subsequently  used  to  constrain 
the  hydrologic  model  in  a  second  inversion.  In  the  coupled  inversion  GPR  traveltimes  were  used 
to  directly  constrain  the  hydrologic  model.  Several  different  cases  were  investigated  to  evaluate 
the  value  of  using  different  GPR  arrivals  as  a  data  constraint,  which  included  (i)  the  ground  wave 
only,  (ii)  reflection  from  the  wetting  front  only,  and  (iii)  both  the  ground  wave  and  reflected 
wave.  A  comparison  of  the  true  hydrologic  parameters  with  those  estimated  by  the  sequential 
versus  coupled  approach  is  given  in  Table  1.  Note  that  in  all  cases  the  traveltime  data  were 
matched  with  a  low  degree  of  misfit. 

Table  1:  Comparison  of  hydrologic  parameter  estimates  obtained  with  sequential  and  coupled 


inversion  given  different  GPR  traveltime  data  as  a  constraint. 


TRUE 

VALUE 

Ground  Wave  +  Reflection 

Ground  Wave  Only 

Sequential 

Coupled 

Sequential 

Coupled 

Oi 

0.05 

0.063 

0.048 

0.071 

0.048 

Os 

0.35 

0.343 

0.349 

0.257 

0.351 

Ks 

0.72 

0.560 

0.785 

0.004 

0.836 

S 

0.98 

0.919 

0.986 

3.573 

0.995 

N 

2.00 

0.985 

4.736 

3.328 

2.336 

For  the  case  where  both  the  ground  wave  and  wetting  front  reflection  are  used  as  data 
constraints,  there  is  a  substantial  amount  of  information  available  to  capture  the  hydrologic 
behavior:  the  ground  wave  constrains  the  water  content  behind  the  wetting  front  while  the 
reflection  provides  independent  information  on  the  velocity  of  the  wetting  front,  which  is 
dependent  on  both  flux  parameters  (K  and  S)  and  storage  potential  (i.e.,  0s-0i).  As  a  result,  both 
the  sequential  and  coupled  inversion  provide  good  estimates  of  the  hydrologic  model  parameters, 
with  the  exception  of  N  (Table  1).  In  contrast,  when  only  the  ground  wave  is  used  as  a 
constraint,  the  independent  information  on  the  wetting  front  movement  is  lost.  The  sequential 
inversion  therefore  yields  large  uncertainty  in  both  the  wetting  front  depth  and  initial  water 
content,  which  translates  into  poor  estimates  of  the  hydrologic  model  parameters  (Table  1).  The 
coupled  inversion,  on  the  other  hand,  continues  to  provide  excellent  estimates  of  the  hydrologic 
model  parameters  (Table  1).  This  is  because  the  hydrologic  model  provides  an  additional 
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constraint  beyond  the  available  geophysieal  data  -  both  the  traveltime  data  and  the  laws  of 
hydrology  must  be  honored  at  the  same  time  in  this  inversion  scheme.  In  this  way  the  hydrologie 
model  acts  as  a  physically-based  regularization  to  stabilize  the  inversion  and  provide  aecurate 
parameter  estimates. 

Conclusions:  This  part  of  the  study  illustrated  that  the  coupled  inversion  seheme  ean 
signifieantly  outperform  sequential  inversion  for  estimating  unsaturated  flow  parameters. 


B.  Laboratory  Investigation  of  GPR  Response  to  Infiltration 


Observed  Data 


Experiment  Time  [min] 

Mean  Water  Content  Between  4 -7 cm  Depth 


Experiment  Time  [min] 


The  second  part  of  this  study 
evaluated  the  response  of  GPR  to  infiltration 
events  using  empirieal  data.  To  this  end  an 
irrigation  system  was  developed  that  eould 
provide  a  speeified  flux  of  water  to  the  upper 
boundary  of  a  sand  box  while  eoineidently 
monitoring  the  infiltration  event  with  GPR. 
The  antennas  were  kept  at  a  fixed  location, 
but  sampled  eontinuously  during  the 
experiment.  In  the  experiments,  the  flux 
applied  to  the  upper  boundary  was  varied  as 
shown  in  Figure  Id  to  ereate  multiple  periods 
of  transient  and  pseudo-steady  state 
conditions.  Water  content  sensors  were 
located  in  the  upper  portion  of  the  box 
(~5.5cm  depth)  and  lower  portion  of  the  box 
(~45cm)  to  independently  monitor  moisture 
changes  during  the  experiment. 

The  GPR  data  obtained  using 
900MHz  antennas  are  shown  in  Figure  la. 
Arrivals  related  to  the  ground  wave  and 
reflections  originating  from  the  bottom  of  the 
tank  are  elearly  seen  in  the  data.  Also 
apparent  are  reflections  related  to  wetting 
fronts  resulting  from  both  the  initiation  of 
flow  at  early  times  and  from  the  first  ehange 
in  flux  after  43  minutes  from  the  start  of  the 
experiment.  Subsequent  wetting  fronts  at 


Experiment  Time  [mm] 


Figure  1:  (a)  Observed  GPR  response  during 
infiltration  event,  dotted  line  shows  comparison  to 
arrival  times  in  simulated  data;  (b)  simulated  GPR 
response  based  on  lab  measurements  of  hydraulic 
parameters,  (c)  average  observed  and  simulated 
water  content  between  4-7cm  depth,  (d)  water 
application  schedule. 
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later  times  appear  to  have  low  amplitude  given  that  the  soil  is  already  significantly  wetted.  The 
patterns  observed  in  the  data  are  in  good  agreement  with  the  water  content  trends  observed  using 
the  soil  moisture  probes  (Figure  Ic).  Similar  results  were  observed  using  450MHz  antennas 
except  that  the  separation  between  the  ground  wave  and  reflections  from  the  wetting  front 
occurred  at  later  times  due  to  the  lower  frequency.  As  a  result,  it  was  more  difficult  to 
discriminate  the  wetting  front  reflections  in  that  data. 

To  evaluate  the  empirical  data,  simulations  of  the  GPR  response  were  performed  using 
the  flow  model  HYDRUS-ID  [Simunek  et  al.,  2008]  and  a  MATLAB  based  2D  GPR  simulator 
created  by  Irving  and  Knight  [2006].  The  flow  was  simulated  using  the  known  boundary 
conditions,  i.e.,  the  applied  flux  in  Figure  Id  at  the  upper  boundary  and  a  seepage  face  condition 
on  the  lower  boundary.  The  Mualem-van  Genuchten  model  was  to  fit  observed  data  for  the 
unsaturated  hydraulic  conductivity  and  pressure-saturation  function.  The  resulting  predicted 
GPR  response  is  shown  in  Figure  lb  and  predicted  water  content  near  the  moisture  probes  in 
Figure  Ic.  The  match  between  the  simulated  GPR  response  and  the  observed  data  is  excellent 
given  that  the  model  prediction  is  based  on  independent  samples  and  not  calibration  to  the 
observed  response.  The  ground  wave,  second  wetting  front,  and  bottom  reflection  are  all  clearly 
visible.  An  exception  is  the  reflection  from  the  first  wetting  front,  which  can  be  found  in  the 
simulated  data  but  has  much  lower  amplitude  than  in  the  empirical  response.  The  traveltime  for 
each  of  the  main  arrivals  in  the  simulated  data  are  plotted  overtop  of  the  empirical  data  in  Figure 
la  to  provide  a  direct  comparison  between  the  results.  It  is  apparent  that  there  is  generally  a 
good  match  except  during  the  period  between  50-75  minutes  where  the  simulated  traveltimes  are 
less  than  the  observed  response.  The  discrepancy  between  the  observed  and  simulated  water 
content  can  be  explained  by  the  underestimation  of  observed  water  content  by  the  simulations  in 
the  near  surface  during  this  period  (Figure  Ic). 

Conclusions:  GPR  is  highly  sensitive  to  infiltration  processes  showing  strong  pattern  responses 
originating  from  the  ground  wave,  reflections  from  wetting  fronts,  and  subsurface 
heterogeneities  (i.e.,  the  tank  bottom  in  this  case).  These  responses  form  trajectories  in  the 
hydrologic  data  that  are  characteristic  of  hydrologic  processes  (e.g.,  soil  wetting/drying). 
Current  simulation  tools  are  sufficient  to  represent  the  GPR  response  to  infiltration  under 
uniform  flow  in  homogeneous  soils. 

C.  Coherency  Analysis  of  GPR  Data  for  the  Estimation  of  Hydrologic  Parameters 

A  challenge  that  was  identified  in  this  research  is  that  the  MCMC  approach  to  estimating 
hydrologic  parameters  from  GPR  data  discussed  above  will  likely  be  difficult  to  generalize  for 
arbitrary  hydrologic  conditions.  This  is  because  models  with  a  sharp  wetting  front,  such  as  the 
Philip  [1957]  model,  are  limited  to  fairly  simple  hydrologic  scenarios  such  as  homogeneous  soils 
with  uniform  initial  water  content  and  simple  boundary  conditions  (i.e.,  fixed  head  or  flux). 
Though  in  some  cases  these  models  are  appropriate  and  the  approach  discussed  in  Section  A  can 
be  applied,  we  found  that  simplified  models  are  not  always  able  to  provide  reliable  simulations 
of  water  content  for  variable  flux  boundaries,  such  as  that  shown  in  Figure  1,  given  arbitrary  soil 
parameters.  As  a  result,  numerical  models  for  flow  must  be  used  to  simulate  water  content  under 
more  general  conditions.  In  this  case,  a  sharp  wetting  front  may  not  be  defined  and  ray  tracing  or 
full  numerical  solutions  to  MaxwelFs  equations  need  to  be  used  to  simulate  GPR  responses. 
The  increased  computational  burden  could  become  a  limiting  factor  for  the  computationally 
intensive  MCMC  techniques.  To  address  this  problem,  a  new  approach  for  estimating 
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hydrologic  properties  that  makes  use  of  the  trajeetories  observed  in  Figure  1  was  hypothesized 
and  investigated  as  a  eomputationally  efficient  means  for  determining  unsaturated  flow 
parameters  under  arbitrary  conditions.  The  approaeh  is  conceptually  similar  to  normal  moveout 
methods  used  to  estimate  wave  veloeity  [e.g.,  Fisher  et  al.,  1992], 

A  common  method  for  analyzing  multi-offset  GPR  data  is  to  calculate  a  measure  of 
coherency  along  the  hyperbolic  trajectory,  i.e.,  normal  moveout,  deseribing  the  ehange  of 
traveltime  with  antenna  offset  as  a  function  of  EM  wave  velocity  [Neidell  and  Taner,  1971]. 
This  method  works  beeause  the  shape  of  the  normal  moveout  trajectory  is  directly  related  to  the 
wave  velocity.  The  veloeity  that  ereates  a  trajectory  through  the  data  along  which  the  similarity 
between  the  traees  is  maximized  (i.e.,  traces  constructively  interfere  when  the  normal  moveout 
effeet  is  taken  into  aeeount)  is  likely  the  true  veloeity. 

In  the  same  way,  it  was  speeulated  in  this  research  that  GPR  data  can  be  analyzed  by 
calculating  coherence  measures  along  hydrologically  defined  trajectories  (e.g..  Figure  1)  to 
determine  hydrologie  parameters  eontrolling  flow.  In  this  work  semblance  and  the  signal  power 
for  data  windows  following  the  trajeetories  were  used  as  measures  of  coherency.  A  key 
problem,  however,  is  to  define  the  hydrologic  trajectories  in  a  computationally  efficient  manner. 
For  situations  where  a  model  with  a  sharp  wetting  front  ean  be  used  this  is  easily  aeeomplished 
using  ray-based  ealeulations.  However,  for  arbitrary  hydrologie  scenarios  an  alternative  method 
is  proposed  for  finding  the  trajeetories  using  GPR  reflection  coefficients.  First,  water  eontents 
are  calculated  with  a  numerical  simulator  (i.e.,  HYDRUS-ID)  as  a  funetion  of  time  and  depth 
during  the  infiltration  event,  whieh  can  then  be  transformed  to  transient  dielectrie  constant 
profiles.  Converting  the  dielectrie  constants  to  velocity  allows  for  the  reflection  coefficients  to 
be  mapped  as  a  funetion  of  GPR  arrival  time.  The  reflection  coefficients  are  then  filtered  to 
identify  diserete  arrivals  that  form  hydrologie  trajeetories  that  are  dependent  on  the  initial 
unsaturated  flow  parameters  (and  boundary  conditions)  used  to  drive  the  infiltration  model.  The 
coherency  of  GPR  data  along  the  trajectories  ean  then  be  evaluated  by  extracting  the  appropriate 
data  from  the  observed  GPR  response.  Therefore,  the  underlying  hypothesis  of  the  approaeh  is 
that  the  true  hydrologie  parameters  will  produee  trajeetories  that  maximize  the  eoherence  of  the 
GPR  data. 


a[x10'^  1/cm] 


^sat 


Figure  2:  Signal  power  along  the  hydrologic  trajectories  as  a  function  of  the  parameters  for  the  Mualem-van 
Genuchten  model.  The  dashed  vertical  line  indicates  the  value  of  the  true  parameter.  The  maximum  signal 
coherence  is  captured  near  the  true  parameter  value  except  for  the  n  parameter  which  does  not  appear  to  have  a 
unique  maximum. 


An  initial  step  toward  testing  this  hypothesis  was  completed  by  performing  a  sensitivity 
analysis  using  the  data  shown  in  Figure  lb.  Figure  2  shows  the  total  signal  power  in  the  GPR 
data  as  a  function  of  each  of  the  hydrologic  parameters.  The  result  suggests  that  coherency  of 
the  signal  is  maximized  by  following  a  trajectory  in  the  GPR  data  that  is  based  on  the  true 
unsaturated  flow  parameters.  The  single  major  exception  is  the  n  parameter  of  the  van 
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Genuchten  model  for  the  water  retention  curve,  which  is  consistent  with  the  result  found  in 
Section  A  where  the  coupled  inversion  of  traveltime  data  provided  a  poor  constraint  on  the  shape 
parameter  (N)  of  the  unsaturated  hydraulic  conductivity  function. 

Conclusions:  Analysis  of  signal  coherency  along  hydrologic  trajectories  in  GPR  data  could 
provide  a  means  to  optimize  hydrologic  model  parameters  with  minimal  computational  effort. 

D.  Field  Observations  of  GPR  Response  to  Infiltration 

Field-scale  experiments  were  conducted  to  perform  an  exploratory  investigation  of 
transient  GPR  responses  during  infiltration  in  heterogeneous  soils.  In  these  experiments,  we 
used  a  sprinkler  system  to  irrigate  an  approximately  10m  x  5m  region  of  a  silty-sand  soil.  GPR 
was  profded  across  the  site  at  multiple  times  during  the  infiltration  experiment.  The  GPR  was 
profiled  between  fixed  stations,  but  due  to  an  inability  to  integrate  real-time  positioning 
information  with  our  current  GPR  system  and  the  discontinuation  of  odometer  wheel  triggering 
devices  by  the  manufacturer  of  our  radar,  we  were  unable  to  precisely  position  the  data.  Instead, 
we  attempted  to  consistently  move  the  antennas  at  a  constant  velocity  across  the  site  - 
introducing  the  potential  for  positioning  errors.  Despite  the  possible  positioning  errors.  Figure  3 
illustrates  significant  shifts  in  the  arrival  time  of  the  ground  wave  and  different  reflections  across 
the  site  as  a  result  of  the  infiltration  event.  In  some  areas,  there  is  a  change  in  the  character  of 
the  GPR  image,  but  it  is  currently  unclear  whether  this  is  an  effect  resulting  from  the  water 
infiltration  or  errors  in  positioning  the  antennas. 


Figure  3:  Two  10m  transects  of  GPR  data  collected  during 
a  field  infiltration  experiment.  The  upper  image  shows  the 
data  collected  prior  to  the  start  of  the  irrigation  sprinklers. 
The  bottom  image  shows  the  same  transect  profiled  after  16 
minutes  of  water  application  to  the  site.  Note  the  shift  in 
fhe  groundwave  and  other  reflections  is  a  result  of 
increased  water  content  (decreased  wave  velocity)  of  the 
soils. 


Conclusions:  Shifts  in  GPR  responses  are  likely  to  be  clearly  observable  in  field  data.  In  this 
study  we  found  that  improved  instrumentation  for  the  positioning  of  the  GPR  antennas  is  needed 
to  improve  the  data  quality  to  allow  for  calibration 
of  hydrologic  models  with  GPR  in  the  field. 
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PROJECT  CONCLUSIONS 


This  short  term  research  project  has  demonstrated  the  potential  of  GPR  for  constraining 
the  parameters  of  unsaturated  flow  models  using  the  coupled  approach  to  hydrogeophysical 
estimation.  This  is  particularly  true  when  GPR  data  contain  insufficient  information  to 
independently  constrain  the  hydrologic  model.  In  this  case,  the  coupled  approach  allows  the 
hydrologic  model  to  regularize  the  inversion  and  produce  reliable  estimates  of  hydrologic 
parameters.  Stochastic  calibration  of  models  using  traveltime  data  was  shown  to  be  effective  for 
simple  flow  scenarios.  Coherency  analysis  along  hydrologic  trajectories,  a  new  methodology 
that  was  developed  as  a  result  of  this  research,  shows  promise  for  identifying  hydrologic  model 
parameters  for  more  complicated  flow  scenarios.  The  behavior  of  the  ground  wave  during 
infiltration  was  one  particular  point  of  interest  in  this  research.  While  it  was  shown  that  the 
ground  wave  response  is  sensitive  to  flow  processes,  this  research  suggests  that  reflection  data 
provide  additional  information  on  the  rate  of  flow  and  can  be  readily  included  in  model 
calibration  algorithms.  Therefore,  it  is  recommended  that  all  GPR  arrivals  be  used  to  calibrate 
hydrologic  models  when  possible. 
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APPENDIX 


STIR:  CHARACTERIZATION  OF  NEAR-SURFACE  MOISTURE  DYNAMICS  USING 
HYDROLOGIC  MODEL-BASED  INTERPRETATION  OF  GPR  DATA 

ARO  Project  Number:  54969-EV-II 

-  FINAL  TECHNICAL  REPORT  - 

Abstract:  This  report  provides  a  brief  summary  of  activities  and  accomplishments  completed 
under  ARO  Project  Number  54969-EV-IL  The  main  objective  of  this  short-term  research  project 
was  to  test  the  hypothesis  that  the  GPR  ground  wave  provides  useful  information  for 
constraining  hydrologic  processes  in  the  vadose  zone.  To  investigate  this  problem  a  series  of 
modeling  and  empirical  studies  were  performed.  The  modeling  studies  examined  whether  GPR 
data  could  be  directly  used  to  constrain  unsaturated  hydraulic  parameters  using  different 
approaches  to  parameter  estimation.  The  coupled  inversion  approach,  where  hydrologic  and 
geophysical  models  are  linked  together  into  a  single  forward  model,  was  found  to  be  an  effective 
strategy  for  estimation  of  the  flow  parameters.  Empirical  studies  focus  on  both  laboratory  and 
field  tests.  Lab  studies  were  performed  in  a  sand  tank  with  spatially  uniform  applied  infiltration. 
The  results  of  these  tests  showed  excellent  sensitivity  of  GPR  to  infiltration  processes  and 
support  the  conclusions  from  the  modeling  studies  that  GPR  can  constrain  hydrologic  processes. 
In  practice,  however,  multiple  arrivals  should  be  used  in  interpretation  since  the  ground  wave 
arrival  is  obscured  by  the  reflection  caused  by  the  wetting  front  at  early  times  during  infiltration 
events.  Field  tests  were  performed  to  determine  if  transient  patterns  in  GPR  data  could  be 
observed  in  a  heterogeneous  environment  under  non-uniform  flow  conditions.  These  tests 
showed  promising  evidence  to  qualitatively  suggest  that  the  methodology  used  in  this  work  will 
be  applicable  to  field  operations.  A  new  method  for  the  analysis  based  on  measuring  coherency 
along  hydrologic  trajectories  in  GPR  data  was  developed  and  shown  to  provide  good  sensitivity 
to  the  parameters  of  the  Mualem-van  Genuchten  model  of  unsaturated  soils. 

1.  Introduction 

The  hypothesis  tested  in  this  research  is  that  hydrologic  model-based  interpretation  of 
transient  GPR  (ground-penetrating  radar)  data,  the  ground  wave  in  particular,  can  lead  to 
improved  characterization  of  soil  water  dynamics.  The  strong  dependence  of  the  high-frequency 
dielectric  permittivity  of  a  soil  on  water  content  has  made  electromagnetic  measurements  (e.g., 
time  domain  reflectometry  (TDR),  ground-penetrating  radar  (GPR),  and  satellite  radars)  popular 
tools  for  monitoring  near  surface  moisture  conditions  [e.g.,  Huisman  et  al.,  2001;  Entekhabi  et 
al.,  2004].  The  GPR  ground  wave  is  of  particular  interest  for  providing  high-resolution  maps  of 
near-surface  water  content  at  the  catchment  scale  [Grote  et  al.,  2003;  Weihermuller  et  al.,  2007]; 
such  data  would  be  invaluable  for  applications  including  understanding  local  rainfall-runoff  and 
infdtration  processes  as  well  as  calibrating  satellite-based  radar  measurements  for  real-time 
monitoring  of  regional  moisture  conditions.  A  limitation  of  current  interpretation  techniques, 
however,  is  the  assumption  that  the  effective  permittivity  derived  from  a  GPR  measurement  can 
be  used  to  estimate  the  average  water  content  of  a  soil,  which  can  subsequently  be  used  to 
interpret  hydrologic  processes.  Previous  studies  have  shown  that  problems  like  preferential 
sampling  of  high  velocity  zones  by  EM  waves  can  lead  to  systematic  biases  in  water  content 
estimates  [Moysey  and  Knight,  2004].  This  problem  is  compounded  if  the  biased  water  contents 
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are  then  used  to  ealibrate  non-linear  unsaturated  flow  models.  Following  this  sequential 
approach  to  geophysical  data  integration,  i.e.,  estimating  transient  water  contents  with  GPR  and 
subsequently  performing  hydrologic  analysis,  may  therefore  lead  to  both  inaccurate  estimates  of 
water  content  and  poor  predictions  of  unsaturated  flow  and  transport. 

This  study  investigated  a  new  methodology  for  characterizing  the  evolution  of  soil 
moisture  using  ground  penetrating  radar.  A  unique  aspect  of  this  study  is  that  the  GPR  data  (e.g., 
travel  times)  were  directly  linked  to  soil  properties  controlling  infiltration  (e.g.,  hydraulic 
conductivity)  by  coupling  hydrologic  process  models  with  geophysical  instrument  models.  We 
demonstrate  this  approach  using  numerical  experiments.  The  sensitivity  of  GPR  to  transient 
infiltration  processes  was  investigated  empirically  using  ID  infiltration  experiments  performed 
in  a  laboratory  sand  tank.  These  experiments  were  used  to  verify  that  specific  arrivals  could  be 
identified  in  the  GPR  data  and  to  evaluate  the  sensitivity  of  the  data  to  infiltration.  Finally,  we 
describe  exploratory  field  experiments  that  were  designed  to  investigate  whether  the  proposed 
methodology  could  be  applied  to  field  studies. 


2.  Background 


airwave 
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2.1  Estimation  of  water  content  using  GPR 

GPR  data  -  the  GPR  ground  wave  in  particular  -  has  received  increasing  attention  as  a 
tool  for  mapping  near  surface  water  content  variations  [Du  and  Rummel,  1994;  Du,  1996;  van 
Overmeeren  et  ah,  1997;  Weiler  et  ah,  1998;  Grote  et  ah,  2003;  Huisman  et  ah,  2001,  2003; 
Klysz  and  Balayssac,  2007;  Strobbia  and  Cassiani,  2007;  Weihermuller  et  ah,  2007].  The 
ground  wave  is  the  portion  of  energy  emitted  from  a  GPR  transmitter  antenna  that  travels 
directly  to  the  receiver  antenna  through  the  soil  (Figure  1).  The  travel  time  between  the 
transmitter  and  receiver  can  be  identified  in  the  radar  data  and  used  to  determine  the  velocity  of 
the  EM  wave,  which  is  primarily  related  to  the  dielectric  permittivity  of  the  soil.  In  most  soils, 

dielectric  permittivity  is  strongly  dependent  on  water 
content  [e.g.,  Topp  et  ah,  1980;  Huisman  et  ah, 
2001].  By  keeping  the  transmitter  and  receiver 
antennas  separated  by  a  fixed  spacing  as  they  are 
moved  across  the  ground  surface  it  is  possible  to 
quickly  map  near-surface  water  contents  over  large 
regions  [e.g.,  Grote  et  ah,  2003]. 

In  addition  to  the  ground  wave,  however, 
multiple  other  arrivals  are  also  produced  as  a  result 
of  alternate  pathways  transferring  energy  between 
source  and  receiver  antennas  (Figure  2).  Strobbia 
and  Cassiani  [2007]  recently  pointed  out  how 
refracted  waves  could  also  severely  impact  the 
interpretation  of  water  content  using  the  GPR 
ground  wave.  Refracted  waves  occur  when  energy 

Near-surface  water  content  ^  P^^h  that  bypasses  low-velocity  zones, 

can  be  identified  from  thereby  potentially  arriving  at  the  receiver  earlier 

fluctuations  in  the  arrival  of  the  GPR  than  direct  or  reflected  waves.  Soil  velocities  can  be 

ground  wave.  The  constant-offset  significantly  overestimated  if  the  arrival  of  a 

measurement  geometry  is  shown  above  refracted  wave  obscures  the  ground  wave  or  is 

example  data  (after  Grote  et  ah,  2003). 


Figure  1: 

variations 
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Figure  2:  (a)  Possible  ray  paths  for 
energy  traveling  between  GPR  transmitter 
and  receiver  antennas  when  a  wet  layer  is 
present  near  the  surface,  (b)  Expected 
arrival  time  for  each  ray  path  as  a 
function  of  antenna  separation.  Note  that 
the  ground  wave  arrival  may  be  obscured 
at  small  offsets  by  the  air  wave  and  at 
larger  offsets  by  refracted  waves,  (c)  An 
example  of  CMP  data  that  can  be  used  to 
optimize  the  transmitter-receiver  survey 
to  minimize  wave  interference  for 
profiling  studies.  (From  Strobbia  and 
Cassiani,  2007). 


^2^  CriticaHy  refracted  wave 


b)  Transmitter-receiver  offset  (m)  c)  Transmitter-receiver  offset  (m) 


mistakenly  identified  as  the  ground  wave.  Given  that  infiltration  events  by  nature  generate  slow 
(i.e.,  wet)  zones  overlying  fast  (i.e.,  dry)  zones,  refracted  waves  might  be  expected  to  be  a 
common  problem.  This  is  a  particularly  serious  issue  for  estimating  water  content  from  fixed 
offset  profiles  (of  the  type  shown  in  Figure  1).  When  multi-offset  profiles  are  collected, 
distinctive  patterns  in  the  data  (e.g.,  Figure  2b)  can  be  used  to  discriminate  between  different 
arrivals  and  thereby  guide  interpretation.  However,  van  der  Kruk  [2006]  recently  discussed  how 
low  velocity  layers,  like  infiltration  fronts,  can  act  as  a  dispersive  wave  guide  that  lead  to  a 
“shingling”  effect  in  multi-offset  GPR  data;  the  implication  being  that  incorrect  velocities  will  be 
estimated  even  in  this  approach  without  special  inverse  procedures.  These  issues  are  significant 
problems  for  fixed  offset  GPR  data.  Discrimination  based  on  the  data  alone  is  not  possible  and 
errors  in  the  estimation  of  the  ground  wave  velocity  are  likely  to  occur.  In  contrast, 
discrimination  between  alternate  subsurface  models  can  be  effectively  achieved  when  the  field 
data  are  compared  to  the  results  of  numerical  simulations,  even  when  interference  from  refracted 
waves  occurs  [e.g.,  Strobbia  and  Cassiani,  2007]  or,  we  hypothesize,  when  the  wave  guide  effect 
is  present. 

Another  significant  challenge  in  using  the  ground  wave  to  investigate  soil  water 
dynamics  is  that  the  support  volume  of  the  measurement  is  not  well  defined.  Grote  et  al.  [2003] 
cite  half  the  Fresnel  zone  as  the  potential  depth  of  investigation,  thus  making  the  sample  volume 
dependent  on  measurement  frequency,  wave  velocity  in  the  soil,  and  antenna  separation.  In  a 
numerical  study  Galagedara  et  al.  [2005b]  specifically  investigated  ground  wave  velocity  in  a 
material  with  a  near  surface  wet  or  dry  layer,  i.e.,  scenarios  typical  of  wetting  and  drying  soils. 
They  found  that  the  depth  of  investigation  was  approximately  60%  of  the  GPR  wavelength,  thus 
confirming  the  dependence  on  measurement  frequency  and  soil  velocity.  More  importantly, 
when  the  layer  thickness  was  between  ~0.25-0.6  of  the  GPR  wavelength,  they  observed  a  nearly 
linear  dependence  of  soil  velocity  on  layer  thickness.  It  is  clear  that  with  this  type  of  simple 
linear  mixing,  i.e.,  Vsoii  =  f  Vi  +  (1-f)  V2,  it  is  straight  forward  to  determine  the  apparent  velocity 
of  a  soil  (vsoii)  if  the  depth  of  a  wetting  front  (f)  (relative  to  the  GPR  wavelength)  is  known  along 
with  the  velocity  (i.e.,  water  content)  within  the  wetted  (vi)  and  dry  (V2)  zones.  However,  it  is 
impossible  to  uniquely  estimate  the  water  content  in  the  wetted  zone  without  making 
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Figure  3:  Illustration  of  the  non-unique  nature  of 
ground  wave  data  for  estimating  the  velocity  and 
thickness  of  a  near-surface  layer  given  a  simple 
linear  mixing  law.  Every  combination  of 
parameters  shown  will  produce  an  effective  (i.e., 
ground  wave)  velocity  of  exactly  0.122  m/ns.  Each 
line  represents  a  different  possible  velocity  for  the 
region  below  the  layer.  Essentially  any  water 
content  could  be  predicted  in  the  near  surface  layer 
unless  additional  constraints  on  water  content 
below  the  layer  and  layer  thickness  are  given. 


assumptions  regarding  the  depth  of  the  wetting 
front  and  initial  water  content  (Figure  3).  Thus 
the  measurements  of  ground  wave  velocity  may 
not  lead  to  a  good  representation  of  moisture 
conditions  in  the  very  near  surface  environment 
(e.g.,  <10cm),  which  is  the  critical  interface 
between  the  subsurface,  land  surface  and 
atmosphere.  This  problem  is  evident  in  the 
empirical  results  of  Galagedara  et  al.  [2005a] 
where  they  found  that  a  good  match  between 
GPR  and  TDR  measurements  during  infiltration 
and  drying  events  could  only  be  obtained  when 
the  length  of  the  TDR  probes  was  similar  to  the 
GPR  investigation  depth. 

The  unknown  support  volume  of  the 
ground  wave,  interference  from  refractions,  and 
preferential  sampling  of  high  velocity  zones  are 
all  problems  that  can  confound  the  estimation  of 
water  content  from  ground  wave  data.  However, 
the  effect  of  each  of  these  issues  can  be  readily 
accounted  for  using  forward  models  describing 
GPR  wave  propagation.  As  a  result,  the 
hydrologic  interpretation  of  ground  wave  data 
could  be  improved  significantly  if  the  analysis  of  the  ground  wave  was  performed  by  comparing 
observed  transient  field  data  with  numerical  simulations  through  inversion.  In  practice,  however, 
even  with  the  use  of  forward  modeling,  a  priori  information  may  be  needed  to  provide  an 
additional  constraint  that  would  allow  investigation  of  subsurface  dynamics  during  infiltration 
with  GPR  (e.g.,  as  demonstrated  by  Figure  2).  This  research  investigates  whether  this  additional 
level  of  constraint  provided  by  transient  data  can  be  imposed  through  the  use  of  hydrologic 
models  to  generate  physically  possible  subsurface  realizations  that  can  be  used  to  drive  GPR 
forward  simulations. 


Relative  Layer  Thickness  [-] 


2.2  Hydrologic  Model-Based  Interpretation  of  Geophysical  Data 

The  current  state-of-the-art  approach  to  integrating  geophysical  methods  into  hydrologic 
problems  is  a  two-step  process  (Figure  4).  First  geophysical  data  are  collected  and  analyzed  to 
extract  information  about  the  distribution  of  geophysical  state  variables  in  the  subsurface  (e.g., 
extraction  of  dielectric  permittivity  from  ground  wave  arrival  times).  Rock  physics  relationships 
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Figure  4:  The  current  state-of-the-art  integration  of  geophysical  data  into  a  hydrologic  problem  follows  a 
sequential  process  where  geophysical  data  are  collected  and  analyzed  to  yield  indirect  estimates  of  hydrologic 
state  variables  (e.g.,  water  content).  In  a  subsequent  analysis  effort,  these  estimates  may  be  used  as  data  within 
an  optimization  scheme  to  calibrate  the  parameters  of  hydrologic  models. 
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are  used  to  eonvert  these  values  to  hydrologic  state  variables  (e.g.,  water  content),  which  may 
subsequently  be  used  to  constrain  the  parameters  governing  hydrologic  process  models. 

A  key  limitation  of  this  approach  is  that  the  resolution  of  geophysical  data  is  often  only 
sufficient  to  uniquely  identify  averages  of  subsurface  properties  (as  demonstrated  in  Figure  3). 
The  inverse  problem  is  therefore  typically  mathematically  stabilized  (i.e.,  regularized)  using  a 
priori  information,  such  as  constraints  on  spatial  continuity  of  the  geophysical  properties  [e.g., 
Menke,  1984].  Often  times  the  choice  of  prior  information  is  arbitrary,  e.g.,  through  the 
application  of  a  smoothness  constraint  on  the  estimated  model,  rather  than  selected  objectively 
based  on  observed  data  or  readily  identifiable  laws.  As  a  consequence,  the  rock  physics 
conversion  from  geophysical  to  hydrologic  properties  at  the  field  scale  can  be  complicated, 
spatially  variable,  and  dependent  on  both  measurement  geometry  and  subsurface  properties 
[Moysey  et  al.,  2005;  Singha  and  Moysey,  2006;  Singha  et  al.,  2007].  Moysey  et  al.  [2006] 
speculated  that  this  non-uniqueness  problem  (and  therefore  the  associated  field-scale  rock 
physics  problem)  could  be  overcome  if  the  laws  of  hydrology  are  used  to  enforce  a  priori 
constraints  on  the  distribution  of  geophysical  state  variables  rather  than  making  arbitrary  choices 
about  the  spatial  distribution  of  these  properties. 

Using  hydrologic  models  to  constrain  the  inversion  of  geophysical  data  (Figure  5)  is  a 
relatively  new  idea  that  has  not  yet  received  an  extensive  amount  of  attention  [Rucker  and  Ferre, 
2004;  Kowalski  et  al.,  2005;  Lambot  et  al.,  2006;  Sicilia  and  Moysey,  2007;  Fowler  and  Moysey, 
2007;  Moysey  et  al.,  2007;  Hinnell  et  al,  2008  {submitted  to  WRR)].  The  approach  initiates  by 
selecting  a  set  of  test  parameters  to  run  a  hydrologic  model  used  to  predict  the  spatial  and 
temporal  distribution  of  hydrologic  state  variables  through  time.  Point-scale  rock  physics 
relationships,  which  are  easier  to  define  than  field-scale  rock  physics  relationships  [Moysey  et 
al.,  2006],  are  used  to  transform  the  hydrologic  state  variables  to  geophysical  properties.  A 
geophysical  instrument  model  is  then  used  to  replicate  the  actual  experiment  performed  in  the 
field.  A  poor  comparison  between  the  simulated  and  observed  field  data  indicates  a  poor  initial 
choice  of  parameters  in  the  hydrologic  model  thereby  leading  to  an  update  of  these  parameters. 
The  model  parameters  are  updated  continuously  until  a  good  match  between  the  simulated  and 
observed  geophysical  data  is  obtained. 
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Figure  5:  Coupling  hydrologic  process  models  with  geophysical  instrument  models  establishes  a  direct  link  between 
geophysical  data  and  hydrologic  parameters.  This  new  approach  to  hydrogeophysical  data  integration  avoids  the 
geophysical  inversion  step  and  therefore  does  not  impose  arbitrary  prior  constraints  on  the  geophysical  properties. 
In  this  approach  a  priori  constraints  on  the  geophysical  properties  are  automatically  imposed  by  the  physical  laws 
built  into  the  hydrologic  model. 


Note  that  in  this  procedure  there  is  no  geophysical  interpretation  step  -  the  geophysical 
properties  are  only  involved  in  forward  models,  not  inverse  models,  so  there  is  no  need  for  the 
assumption  of  a  priori  constraints  on  these  properties.  Also,  there  is  no  restriction  on  the 
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complexity  of  the  geophysieal  data  used  in  the  inversion.  For  example,  rather  than  extraeting 
ground  wave  travel  times  from  the  GPR  data,  the  comparisons  between  simulated  and  observed 
data  could  made  for  the  GPR  amplitude  record.  As  a  result,  interaetions  between  different 
arrivals  may  be  eonsidered  as  an  opportunity  to  identify  and  discriminate  processes,  rather  than  a 
problem  that  reduces  data  quality. 

To  date  few  studies  have  used  the  eoupled  interpretation  approaeh  to  ealibrate  infiltration 
models  using  GPR  data.  Kowalski  et  al.  [2005]  illustrated  that  a  eoupled  analysis  could  be  used 
to  estimate  flow  properties  using  cross-borehole  GPR  monitoring  of  an  infiltration  event. 
Lambot  et  al.  [2006]  performed  a  synthetic  experiment  to  demonstrate  that  air-launehed  radar 
measurements  could  theoretically  be  used  to  constrain  infiltration  processes.  Sicilia  and  Moysey 
[2007]  performed  the  first  eomparative  analysis  of  the  sequential  (Figure  4)  and  eoupled  (Figure 
5)  data  integration  methods.  These  authors  investigated  whether  intrinsie  permeability  values 
could  be  constrained  by  monitoring  an  infiltration  event  with  eross-borehole  GPR.  They  found 
that  the  sequential  analysis  resulted  in  permeability  estimates  that  consistently  underestimated 
true  values  by  an  order  of  magnitude  (Figure  6a).  In  contrast,  the  coupled  analysis  resulted  in  no 
permeability  bias  and  overall  produced  mueh  better  fits  between  simulated  and  observed  data 
(Figure  6b).  Hinnell  et  al.  [2008  (submitted  to  WRR)]  present  a  synthetie  study  where  surfaee- 
based  eleetrical  resistivity  monitoring  of  ID  infiltration  is  used  to  estimate  5  different  soil 
properties.  These  authors  used  an  implementation  of  a  Markov  Chain  Monte  Carlo  method  to 
compare  both  the  estimation  aceuraey  and  uncertainty  resulting  from  the  sequential  versus 
coupled  analysis  strategies.  Figure  7  shows  predietions  of  the  wetting  front  depth  versus  time 
based  on  the  flow  models  ealibrated  with  the  resistivity  measurements.  Note  that  the  data  used 
in  either  ease  is  identical;  the  only  differenee  between  the  results  is  whether  the  sequential  or 
coupled  interpretation  teehnique  was  used.  These  studies  illustrate  the  potential  performanee 
improvement  that  can  be  gained  by  enforeing  hydrologic  constraints  on  the  geophysieal  data. 
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<— Figure  6:  Comparison  of 
objective  functions  for  the 
estimation  of  permeability 
using  (a)  sequential  and  (b) 
coupled  analyses  of  a  GPR 
monitored  infiltration  test. 
Note  that  the  sequential 
approach  results  in  a  bias  in 
the  estimate  and  higher  data 
misfit  than  the  coupled 
approach  (Sicilia  and 
Moysey,  2007). 

^Figure  7:  Predicted 
depth  of  an  infiltration  front 
using  a  model  calibrated 
with  resistivity  data  using 
(a)  sequential  and  (b) 
coupled  analyses.  Note 
that  the  sequential  approach 
yields  much  larger  levels  of 
prediction  uncertainty,  yet 
the  true  front  still  falls 
outside  the  95%  confidence 
bounds  of  the  prediction. 
(Hinnell  et  al.,  2008). 
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3.  Overview  of  Methods  and  Results  from  this  Research 

This  research  investigated  the  use  of  GPR  dat  for  constraining  hydrologic  models  using  three 
approaches:  (1)  numerical  modeling,  (2)  laboratory  studies,  and  (3)  field  studies.  The  methods 
and  results  obtained  from  each  of  these  is  outlined  below  in  the  three  following  sections. 


3.1  Numerical  Modeling  Study 


3.1.1  Description  of  Model 

The  modeling  in  this  study  focused  on  comparing  the  sequential  (Figure  4)  and  coupled 
(Figure  5)  inversion  approaches  for  estimating  hydrologic  properties  from  GPR  data.  The  initial 
set  of  numerical  experiments  considered  a  scenario  of  ID  infiltration  under  zero  ponding  depth 
in  a  uniform  soil  with  constant  initial  water  content  (0i).  In  this  case,  infiltration  can  be 
described  using  the  Philip  [1957]  model  where  cumulative  infiltration  I(t)  at  an  arbitrary  time  t  is 

related  to  the  sorptivity  (S)  and  saturated  hydraulic  conductivity  (Ks)  of  the  soil: 

1 

l{t)  =  SD+K^t  (1) 


The  wetting  front  Zwf  in  the  Philip  model  is  a  sharp  boundary  such  that  the  depth  of  the  interface 
can  be  defined  at  any  time  by  the  storage  capacity  of  the  soil  (i.e.,  0wr0i)  and  cumulative 
infiltration: 
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where  the  water  content  above  the  wetting  front  is  0wf,  which  is  equal  to  the  saturated  water 
content  0s  throughout  the  infiltration  period. 

During  the  redistribution  phase  was  modeled  as  a  zero  flux  boundary  and  flow  was 
assumed  to  be  gravity  driven.  In  this  case,  further  migration  of  the  wetting  front  must  result 
from  a  decrease  in  storage,  i.e.,  change  of  water  content  above  the  wetting  front.  We  used  a 
rectangular  drainage  model  [e.g..  Jury  and  Horton,  2004]  to  represent  this  process: 


^wf  - 


(3) 


where  Za  is  the  depth  reached  by  the  wetting  front  at  the  end  of  the  infiltration  period  and  ta  is  the 
time  at  which  infiltration  ceased.  Under  gravity  driven  flow  the  flux  q  is  controlled  by  the 
unsaturated  hydraulic  conductivity  K(0),  i.e.,  q  =  K(0).  In  this  study,  an  exponential  relationship 
was  used  for  the  unsaturated  conductivity  function: 
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where  0  is  the  water  content  at  an  arbitrary  time  t,  0s  is  the  water  content  at  saturation,  and  N  is  a 
soil-specific  parameter.  The  water  content  behind  the  wetting  front  at  any  time  after  the 
cessation  of  infiltration,  i.e.,  t-td,  will  then  be  given  by: 
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Though  this  model  is  simple,  it  is  also  powerful.  For  example.  Figure  8  compares  the  soil 
moisture  distribution  predicted  by  this  model  and  the  numerical  solution  of  Richard’s  equation  as 
simulated  by  the  numerical  model  HydruslD.  Though  the  details  are  slightly  different,  there  is 
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Figure  8:  Comparison  of  the  evolution  of  water  content  profiles  during  an  infiltration  event  as 
determined  by  (a)  solution  of  Richards’  equation  with  Hydrus-ID  and  (b)  the  Philip  infiltration  and 
rectangular  redistribution  model.  Although  the  details  of  the  depth  profiles  are  slightly  different,  the 
simple  model  captures  the  gross  changes  in  water  content  that  are  likely  to  affect  GPR  measurements  in  a 
fraction  of  the  computational  time. 

an  excellent  reproduction  of  the  overall  infiltration  and  redistribution  process.  Furthermore,  the 
model  used  here  can  be  used  to  simulate  flow  in  a  fraction  of  the  time  required  to  run  a  full 
numerical  solution,  thus  is  more  appropriate  for  use  in  stochastic  estimation  algorithms;  the 
simulation  in  Figure  8a  required  18  seconds  of  computational  time  on  a  1.73GHz  desktop 
computer  with  2GB  of  RAM,  whereas  the  simulation  in  Figure  8b  required  only  0.01  seconds  on 
the  same  computer. 

The  results  shown  in  Figure  8  also  emphasize  the  fact  that  the  infiltration-redistribution 
model  conceptualizes  the  subsurface  as  a  time-variable,  single  layer  system,  i.e.,  the  region 
behind  the  wetting  front  above  a  homogeneous  background.  This  can  therefore  also  be  modeled 
as  a  single  layer  system  for  GPR  where  we  have  used  Topp’s  equation  to  relate  water  content 
and  dielectric  constant  [Topp  et  ah,  1980].  Assuming  a  low-loss  environment,  the  velocity  of  the 
EM  wave  is  given  by  v=c/k*’  where  c  is  the  speed  of  light  and  k  is  the  dielectric  constant  of  the 
soil.  The  GPR  travel  time  for  various  pathways  can  be  calculated  analytically  for  this  single 
layer  model  using  ray  theory.  For  example,  the  travel  time  of  the  ground  wave  (Tg)  would 
simply  be  Tg  =  L/vgw,  where  L  is  the  antenna  separation  and  Vwf  is  the  velocity  of  the  EM  wave 
behind  the  wetting  front.  Note  that  at  late  times,  when  the  wetting  front  has  passed  the  sampling 
depth  of  the  ground  wave  Vgw  would  simply  be  the  velocity  of  the  soil  behind  the  wetting  front; 
before  this  time,  however,  the  effective  velocity  will  be  an  average  of  the  velocities  above  and 
below  the  wetting  front  [e.g.,  Galagedara  et  ah,  2005b].  Likewise,  the  travel  time  for  a  reflection 
from  the  wetting  front  would  be  given  by  Twf  =  2*(  (L/2)^  +  Zw/  )'^'^/Vwf. 

3.1.2  Simulation  Details  and  Estimation  Procedure 

To  compare  the  performance  of  the  sequential  and  coupled  inversion  strategies,  a 
numerical  experiment  was  conducted  where  infiltration  was  simulated  for  5  minutes  followed  by 
a  28.5  minute  period  of  redistribution.  The  “true”  soil  properties  used  for  this  test  are  given  in 
Table  1.  The  resulting  water  content  distribution  as  a  function  of  time  is  shown  in  Figure  9.  The 
GPR  travel  times  for  the  ground  wave  and  reflection  from  the  wetting  front  boundary  are  given 
in  Figure  10;  1%  random  Gaussian  noise  was  added  to  each  travel  time. 


Table  1 
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igure  9: 


Position  of  the  wetting  front  and  water  eontents  as  a  function  of 
time  during  the  simulated  infiltration  experiment. 


Figure  10:  Travel  times  for  the  GPR  ground  wave  and  reflection  from  the 
wetting  front  during  infiltration  and  redistribution.  Random 
Gaussian  noise  (1%)  was  added  to  each  travel  time. 
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The  noisy  travel  times  in  Figure  10  represent  the  data  used  to  eonstrain  the  hydrologie 
parameters  of  the  infiltration  and  redistribution  model  using  both  the  sequential  and  eoupled 
inversion.  In  the  case  of  sequential  inversion  (Figure  4),  the  travel  times  were  first  used  to 
estimate  EM  wave  velocity  above  of  the  wetting  front  and  the  depth  to  the  interface.  The 
velocities  were  then  converted  to  water  contents  and  used  as  a  constraint  in  a  second 
optimization  to  determine  the  hydrologic  properties.  For  the  coupled  inversion,  in  contrast,  the 
travel  times  were  used  as  a  direct  constraint  on  the  hydrologic  model  using  the  approach  outlined 
in  Figure  5.  In  both  cases,  the  Metropolis-Hastings  algorithm  was  used  to  sample  the  posterior 
probability  distribution  of  the  model  parameters.  Gaussian  distributions  were  assumed  for  the 
likelihood  function  in  all  cases  and  candidate  distributions  for  parameters  were  limited  to  the 
ranges  given  in  Table  2. 


Table  2:  Bounds  and  initial  value  of  parameters  used  in  Metropolis-Hastings  algorithm 


Hydrologic  Parameters 

Geophysical  Parameters 

e,  [-1 

0s  [-1 

Ks  [m/day] 

S  [m/day^’^] 

N[-] 

Vwf  [m/ns] 

Vi  [m/ns] 

Zwf[m] 

Lower  Bound 

0 

0.1 

0 

0 

0 

0 

0.08 

0 

Upper  Bound 

0.2 

0.6 

5 

5 

5 

0.1 

0.18 

1 

Initial  Value 

0.01 

0.3 

1.5 

0.5 

1.5 

0.05 

0.18 

0.5 

3.1.3  Summary  of  Results 

Using  MCMC  to  sample  the  posterior  distribution  produces  a  collection  of  viable 
parameter  sets.  For  the  sequential  inversion  10,000  sets  of  the  three  geophysical  variables 
(velocity  above  and  below  the  wetting  front  and  depth  to  the  wetting  front)  were  generated  with 
the  Metropolis-Hastings  algorithm  at  each  observation  time;  the  mean  and  standard  deviation  for 
each  parameter  are  shown  in  Figure  12.  For  both  the  sequential  and  coupled  inversion,  the 
Metropolis-Hastings  algorithm  was  used  to  sample  20,000  samples  of  the  5  hydrologic 
parameters  (i.e.,  0i,  0s,  Ks,  S,  N);  the  mean  and  standard  deviation  of  these  results  are  given  in 
Tables  3  and  4.  For  the  sequential  inversion  this  was  accomplished  by  converting  each  set  of 
velocities  to  water  content  using  the  low-loss  relationship  between  velocity  and  dielectric 
constant  (k  =  (0.3/v)^,  where  v  is  velocity  in  m/ns)  and  then  applying  the  Topp  equation  [1980]. 
The  resulting  water  contents  were  then  used  as  data  to  constrain  the  hydrologic  model.  In 
contrast,  the  GPR  traveltimes  are  used  directly  as  a  constraint  on  the  hydrologic  parameters  for 
the  linked  hydrologic  and  geophysical  model  in  the  coupled  inversion. 

The  fit  of  the  GPR  traveltimes  using  the  sequential  estimation  approach  is  excellent  as 
shown  in  Figure  11.  Despite  the  excellent  fit  to  the  traveltime  data.  Figure  12  indicates  an 
inability  to  reliably  constrain  the  velocity  below  the  wetting  front.  Regardless,  the  velocity 
above  the  wetting  front  and  the  depth  to  the  front  are  estimated  with  reasonable  accuracy.  This  is 
not  surprising  since  for  a  2  layer  model  the  ground  wave  arrival  provides  a  constraint  on  velocity 
above  the  wetting  front  and  the  reflected  wave  provides  an  independent  constraint  on  the  depth 
to  the  interface.  However,  there  is  also  the  potential  for  tradeoffs  between  layer  thickness  and 
apparent  velocity  above  the  wetting  front  as  described  in  Figure  3.  This  effect  is  thought  to  be 
responsible  for  the  “spikes”  apparent  in  the  velocity  above  the  wetting  front  which  are  paired  to 
similar  spikes  in  wetting  front  depth  -  note  that  the  traveltime  data  are  still  fit  very  well  despite 
the  inaccuracy  of  the  velocity  and  depth  estimate. 
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Figure  11:  Fit  of  the  GPR  traveltime  data  using  the  sequential  inversion. 


Figure  12:  Sequential  estimates  of  model  parameters  using  MCMC  at  each  time  of  the  simulation 
-  (a)  velocity  above  and  below  the  wetting  front,  (b)  water  content  above  and  below  the  wetting 
front,  (c)  depth  to  the  wetting  front.  Dots  represent  the  mean  of  all  parameter  sets  accepted  using 
the  Metropolis-Hastings  algorithm  and  dashed  lines  represent  la.  Open  circles  show  the  water 
content  and  wetting  front  depth  predicted  by  the  hydrologic  parameter  using  the  best  fit  (i.e., 
mean)  parameters  from  Table  3.  The  solid  lines  show  the  true  model  response. 
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Table  3:  Estimation  results  using  GPR  traveltimes  from  both  the  ground  wave  and  wetting  front 
reflection.  Mean  and  standard  deviation  of  the  hydrologic  parameters  estimated  from 
samples  from  the  posterior  distributions  using  MCMC.  _ 


TRUE 

VALUE 

SEQUENTIAL  ESTIMATE 

COUPLED  ESTIMATE 

Mean 

a 

Mean 

a 

01 

0.05 

0.063 

0.007 

0.048 

0.001 

Os 

0.35 

0.343 

0.003 

0.349 

0.001 

Ks 

0.72 

0.560 

0.073 

0.785 

0.056 

S 

0.98 

0.919 

0.026 

0.986 

0.007 

N 

2.00 

0.985 

0.378 

4.736 

0.165 

The  fit  to  the  traveltime  data  for  the  coupled  inversion  is  shown  in  Figure  13.  The  main 
difference  from  the  sequential  result  is  that  the  noise  in  the  data  are  not  fit,  while  the  main  trends 
are.  This  is  a  consequence  of  the  fact  that  the  coupled  hydrogeophysical  model  provides  a 
regularizing  effect  on  the  estimation  -  the  noise  in  the  data  cannot  be  explained  by  the 
hydrologic  model,  thus  cannot  be  fit  by  the  coupled  geophysical  and  hydrologic  models.  In  other 
words,  the  coupled  inversion  is  more  stable  and  less  sensitive  to  noise  than  the  sequential 
inversion.  The  hydrologic  response  predicted  by  the  best  fit  parameter  estimates  for  the  coupled 
inversion  are  given  in  Figure  14.  The  match  between  the  true  response  and  that  predicted  by  the 
calibrated  model  is  excellent. 

The  estimation  procedure  described  above  was  repeated  using  only  the  ground  wave 
traveltimes  and  again  using  only  the  traveltimes  for  the  reflection  from  the  wetting  front;  the 
results  are  given  in  Tables  4  and  5,  respectively.  When  only  the  ground  wave  is  used,  the 
velocity  of  the  upper  layer  is  well  constrained,  but  the  constraint  on  the  wetting  front  depth  that 
had  previously  been  provided  by  the  reflection  data  is  lost  (compare  Figure  12c  to  Figure  16c). 
The  results  of  the  sequential  inversion  suggest  that  knowing  the  water  content  of  the  upper  layer 
alone  is  not  sufficient  to  constrain  the  model  parameters.  In  constrast,  in  the  coupled  inversion 
good  estimates  of  the  hydrologic  parameters  can  still  be  achieved  because  the  hydrologic  model 
implicitly  places  an  additional  constraint  of  how  the  water  contents  in  the  upper  layer  are 
allowed  to  change  through  time.  This  extra  information  along  with  the  ground  wave  data 
appears  to  be  sufficient  to  constrain  the  behavior  of  the  model  (Figure  17).  The  contrast  between 
the  performance  of  the  sequential  and  coupled  inversion  schemes  in  this  example  highlights  the 
concept  that  the  hydrologic  model  provides  a  physically  based  regularization  on  the  inversion. 


Table  4:  Estimation  results  using  GPR  traveltimes  from  only  the  ground  wave. 


TRUE 

VALUE 

SEQUENTIAL  ESTIMATE 

COUPLED  ESTIMATE 

Mean 

a 

Mean 

a 

01 

0.05 

0.071 

0.005 

0.048 

0.001 

0s 

0.35 

0.257 

0.001 

0.351 

0.002 

Ks 

0.72 

0.004 

0.004 

0.836 

0.077 

s 

0.98 

3.573 

0.098 

0.995 

0.014 

N 

2.00 

3.328 

0.228 

2.336 

0.589 

Table  5:  Estimation  results  using  GPR  traveltimes  from  only  the  wetting  front  reflection. 


TRUE 

VALUE 

SEQUENTIAL  ESTIMATE 

COUPLED  ESTIMATE 

Mean 

a 

Mean 

a 

01 

0.05 

0.082 

O.OII 

0.055 

0.022 

0s 

0.35 

0.297 

0.015 

0.349 

0.003 

Ks 

0.72 

0.102 

0.093 

0.696 

0.152 

S 

0.98 

0.624 

0.061 

0.964 

0.071 

N 

2.00 

0.3.43 

1.280 

0.822 

0.166 
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Figure  14:  Predicted  response  of  the  hydrologic  model  calibrated  using  coupled  inversion: 

(a)  velocity  above  and  below  the  wetting  front,  (b)  water  content  above  and  below  the  wetting 
front,  (c)  depth  to  the  wetting  front.  Solid  lines  represent  the  true  modeled  response,  whereas  the 
open  circles  show  the  water  content  and  wetting  front  depth  predicted  by  the  calibrated  hydrologic 
model  using  the  best  fit  (i.e.,  mean)  parameters  from  Table  3. 
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Figure  15:  Fit  of  the  ground  wave  traveltime  data  using  the  (a)  sequential  and  (b)  coupled  inversion. 
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Figure  16:  Sequential  estimates  of 
model  parameters  -  (a)  velocity 
above  and  below  the  wetting  front, 
(b)  water  content  above  and  below 
the  wetting  front,  (c)  depth  to  the 
wetting  front. 

Dots  represent  the  mean  of  all 
parameter  sets  accepted  using  the 
Metropolis-Hastings  algorithm  and 
dashed  lines  represent  1  standard 
deviation.  Open  circles  show  the 
water  content  and  wetting  front 
depth  predicted  by  the  hydrologic 
parameter  using  the  best  fit  (i.e., 
mean)  parameters  from  Table  4. 
Solid  lines  show  the  true  model 
response. 


Figure  17:  Predicted  response  of 
the  hydrologic  model  calibrated 
using  coupled  inversion: 

(a)  velocity  above  and  below  the 
wetting  front,  (b)  water  content 
above  and  below  the  wetting  front, 
(c)  depth  to  the  wetting  front. 
Solid  lines  represent  the  true 
modeled  response,  whereas  the 
open  circles  show  the  water 
content  and  wetting  front  depth 
predicted  by  the  calibrated 
hydrologic  model  using  the  best  fit 
(i.e.,  mean)  parameters  from  Table 
4. 
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Figure  18:  Fit  of  the  traveltime  for  the  wetting  front  reflection  using  the  (a)  sequential  and  (b)  coupled  inversion. 


Figure  19:  Sequential  estimates  of 
model  parameters  -  (a)  velocity 
above  and  below  the  wetting  front, 
(b)  water  content  above  and  below 
the  wetting  front,  (c)  depth  to  the 
wetting  front. 

Dots  represent  the  mean  of  all 
parameter  sets  accepted  using  the 
Metropolis-Hastings  algorithm  and 
dashed  lines  represent  1  standard 
deviation.  Open  circles  show  the 
water  content  and  wetting  front 
depth  predicted  by  the  hydrologic 
parameter  using  the  best  fit  (i.e., 
mean)  parameters  from  Table  5. 
Solid  lines  show  the  true  model 
response. 
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Figure  20:  Predicted  response  of 
the  hydrologic  model  calibrated 
using  coupled  inversion: 

(a)  velocity  above  and  below  the 
wetting  front,  (b)  water  content 
above  and  below  the  wetting  front, 
(c)  depth  to  the  wetting  front. 
Solid  lines  represent  the  true 
modeled  response,  whereas  the 
open  circles  show  the  water 
content  and  wetting  front  depth 
predicted  by  the  calibrated 
hydrologic  model  using  the  best  fit 
(i.e.,  mean)  parameters  from  Table 
5. 
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The  results  of  the  sequential  and  eoupled  inversions  using  only  traveltimes  for  the  reflection 
from  the  wetting  front  are  shown  in  Figures  18-20.  In  this  case,  there  is  a  strong  tradeoff 
between  the  velocity  above  the  wetting  front  and  the  depth  to  the  wetting  front  -  many 
combinations  of  velocity  and  depth  can  result  in  the  same  traveltime.  This  fact  is  evident  in  the 
sequential  results,  where  there  is  a  very  high  correlation  between  the  velocity  above  the  wetting 
front  (Figure  19a,  lower  lines)  and  wetting  front  depth  (Figure  19c).  Although  the  sequential 
inversion  converged  to  parameter  reasonable  parameter  estimates  in  this  case,  the  highest  level  of 
uncertainty  was  also  achieved  indicating  the  reduced  degree  of  constraint  of  the  model. 

3.1.4  Conclusions 

The  coupled  inversion  scheme  was  shown  to  provide  better  estimates  of  hydrologic 
model  parameters  compared  to  the  sequential  approach.  Particularly  important  was  the  ability  of 
the  coupled  approach  to  provide  accurate  parameter  estimates  even  under  limited  data  constraint. 
This  result  suggests  that  the  coupled  inversion  approach  regularizes  the  inverse  problem  using 
the  underlying  physics  driving  flow.  Overall,  the  coupled  approach  is  likely  to  provide  better 
results  than  sequential  inversion,  but  due  to  the  regularizing  role  the  hydrologic  model  plays  in 
the  inversion,  further  effort  should  be  made  to  investigate  the  impact  of  conceptual  errors  in  the 
hydrologic  model  on  the  parameter  estimates. 

3.2  Laboratory-Scale  Study  of  GPR  Response  to  Infiltration 

The  second  part  of  the  study  investigated  the  actual  response  of  GPR  to  infiltration 
events.  The  objectives  of  this  work  were  to  evaluate  if  strong  signal  responses  could  be  observed 
from  hydrologically  related  arrivals  (i.e.,  ground  wave  and  wetting  front  reflections).  This 
section  also  describes  a  new  approach  to  the  analysis  of  GPR  data  developed  through  the  course 
of  this  research. 

3.2.1  Experimental  Methods 

The  goal  of  this  study  was  to  generate  a 
uniform  vertical  flow  field  over  the  region  that 
monitored  by  GPR.  The  experiments  were 
performed  in  tanks  filled  with  homogeneous, 
medium  grained  commercial  sand  obtained  from 
CEMEX  USA  (e.g..  Figure  21).  Gravel  was 
placed  under  the  sand  to  allow  for  drainage. 

The  flux  of  water  applied  at  the  top  of  the  tank 
was  controlled  using  a  variable  rate  peristaltic 
pump  (Preston  Monostat,  Cole-Parmer) 
calibrated  for  a  custom  drip  irrigation  system 
with  outlets  created  on  a  1cm  x  1cm  grid  over 
the  entire  tank  surface  by  using  a  syringe  to 
puncture  holes  in  6mm  Tygon  tubing.  Figure  22 
illustrates  the  irrigation  system.  Calibration  experiments  showed  that  spatially  uniform  flow 
could  be  approximately  achieved  using  the  system.  Capacitance  probes  (ECS,  Decagon  Devices, 
Inc.)  were  distributed  at  two  depths  and  at  seven  different  locations  in  the  tank  (shown  in  Figure 
23)  to  monitor  the  uniformity  of  changes  in  near-surface  water  content  and  the  arrival  of  the 
wetting  front  at  the  bottom  of  the  tank.  Tensiometers  were  also  included  to  monitor  pressure  in 


Figure  21;  Example  of  the  packed  sand  boxes  used 
for  the  infiltration  tests. 
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some  experiments,  but  are  not  reported  on  here  as  technical  problems  were  experienced  in 
obtaining  reliable  data. 


The  GPR  system  used  in  this  study  was  a  Sensors 
and  Software  pulseEKKO  1000  with  450  and  900  MHz 
antennas.  The  antennas  were  placed  in  the  center  of  the 
tank  directly  on  top  of  the  tubing  used  for  water 
application.  The  antenna  separation  was  37cm,  as 
measured  from  the  centers  of  the  transmitter  and  receiver 
antenna.  The  antennas  were  not  moved  at  any  point 
during  the  experiments.  GPR  traces  were  collected 
approximately  every  3  seconds  throughout  the 
experiment.  No  processing  was  performed  on  the  GPR 
data. 

Independent  measurements  of  soil  hydraulic 
properties  were  performed  for  samples  of  the  sand  used 
to  pack  the  flow  tank.  The  saturated  hydraulic 
conductivity  (Kg)  of  the  sand  was  measured  using  a 
constant  head  permeameter.  The  saturated  water  content 
(0s)  was  obtained  gravimetrically  by  oven  drying  a 
saturated  sample  of  known  volume.  The  water  retention 
curve  was  determined  using  the  hanging  column  method, 
which  was  fitted  using  the  model  of  van  Genuchten 
[1980]  to  estimate  the  residual  water  content  (0r)  and  parameters  a  and  n.  The  resulting 
parameter  values  from  these  measurements  are  given  in  Table  6. 


Figure  21:  Illustration  of  the  irrigation 
system  with  the  GPR  antennas  resting  on 
the  tubing. 


Figure  23:  Illustration  of  the 
moisture  probe  distribution  used 
to  monitor  uniformity  of  flow 
across  the  sand  tanks.  One  probe 
was  located  in  each  quadrant  of 
the  tank  and  three  probes  were 
located  in  the  center  of  the  tank  to 
identify  non-uniform  flow. 


Table  6:  Measured  hydrologic  properties  of  the  sand  used  in  the  flow  experiments. 


Or 

0s 

* 

a 

n 

Ks 

[cm/min] 

0.041 

0.385 

0.033 

6.26 

1.43 

*a  and  n  are  fitting  parameters  from  the  van  Genuchten 
water  retention  model 
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Initial  experiments  while  the  system  was  being  designed  and  tested  were  performed  using 
the  small  flow  tank  shown  in  Figure  21,  whieh  was  only  30em  deep.  Because  we  dried  the  sand 
between  experiments,  the  use  of  the  smaller  tank  allowed  us  to  cycle  between  wet  and  dry  sand 
during  the  testing  phase  of  the  project.  Figure  24  is  included  here  to  show  the  reason  why  a 
larger  tank  was  used  in  later  experiments  that  we  report  on  in  the  results  section  and  also  to 
illustrate  the  impacts  that  near-surface  heterogeneity  may  have  on  GPR  data.  In  this  data,  the 
tank  depth  is  too  small  to  allow  for  the  separation  of  arrivals  resulting  from  the  movement  of  the 
wetting  front  versus  the  reflection  from  the  bottom  of  the  tank.  Therefore,  it  is  not  possible  to 
uniquely  identify  the  response  of  the  wetting  front.  For  example,  it  is  not  possible  to  pick  the 
traveltime  for  the  reflection  from  the  wetting  front  that  would  be  needed  to  perform  the  type  of 
model  calibration  discussed  in  the  previous  section  of  this  report. 


Start  \vater 


Time  Since  Start  of  Experiment  [min] 


Figure  24:  Example  of  GPR  data  collected  in  a  small  tank  of  30cm  depth.  Note  that  it  is  not 
possible  to  identify  unique  arrivals  related  to  the  wetting  front  in  the  tank  from  those  produced  by 
the  interface  at  the  bottom  of  the  tank.  This  figure  illustrates  why  a  deep  tank  is  required  in  this 
research. 

The  main  experiments  were  conducted  in  a  flow  tank  packed  with  sand  to  a  depth  of 
50cm,  under  which  16cm  of  gravel  was  placed  to  ensure  good  drainage.  During  the  experiments, 
the  rate  of  flow  was  changed  over  time  to  create  both  transient  and  pseudo-steady  state  moisture 
conditions  in  the  tank.  Capacitance  probes  were  distributed  at  depths  of  5.5  and  45  cm  at  the 
same  seven  spatial  locations  in  the  tank  as  shown  in  Figure  23. 

3.2.2  Experimental  Results 

The  results  for  the  flow  experiment  monitored  using  900MHz  and  450MHz  GPR 
antennas  are  shown  in  Figures  25  and  26,  respectively.  The  experiment-specific  schedule  for  the 
applied  flux  to  the  tank  surface  is  given  in  Table  7.  A  variety  of  arrivals  can  be  clearly  identified 
in  the  data,  including  the  air  wave,  ground  wave,  reflections  from  two  different  wetting  fronts, 
and  a  strong  reflection  at  the  boundary  between  sand  and  gravel.  Coherent  noise  at  late  times  is 
likely  the  result  of  reflections  from  the  sides  of  the  tank  and  reflection  multiples. 
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Table  7:  Schedule  of  applied  flux  for  900MHz  experiment;  experiment  ended  after  156.5  minutes. 


Experiment  Time  [min] 

0.0 

4.0 

42.0 

42.3 

64.9 

70.7 

77.5 

77.9 

90.0 

Applied  Flux  [cm/hr] 

0 

2.1 

0 

10 

0 

10 

0 

21 

0 

Experiment  Time  [min] 

95.0 

110.3 

122.7 

131.3 

122.7 

131.3 

136.2 

142.8 

156.5 

Applied  Flux  [cm/hr] 

21 

0 

41 

0 

41 

0 

41 

0 

END 
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Figure  24:  Interpreted  900MHz  GPR  data  during  an  infiltration  experiment. 
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Figure  25:  Interpreted  450MHz  GPR  data  during  an  infiltration  experiment. 
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Table  8:  Schedule  of  applied  flux  for  450MHz  experiment. 


Experiment  Time  [min] 

0.0 

2.0 

40.0 

68.0 

88.0 

95.0 

128.0 

142.0 

Applied  Flux  [cm/hr] 

0.0 

2.4 

10.9 

21.7 

0.0 

38.8 

0.0 

END 

There  are  several  significant  differences  in  the  results  for  the  two  frequencies.  Perhaps 
most  importantly,  the  reflection  from  the  initial  wetting  front  is  obvious  in  the  900MHz  data, 
whereas  it  is  difficult  to  identity  in  the  450MHz  data  because  of  interference  between  it,  the 
ground  wave  and  the  reflection  from  the  bottom  of  the  tank.  Another  difference  is  in  the 
response  of  the  ground  wave  arrival  time.  There  is  a  much  greater  change  in  the  ground  wave 
traveltime  at  early  times  during  the  experiment  for  the  900MHz  data  than  the  450MHz  data. 
This  type  of  behavior  is  consistent  with  the  difference  in  the  sampling  depth  of  the  ground  wave 
for  the  two  frequencies;  at  lower  frequency  the  sampling  depth  is  larger  resulting  in  a  less 
pronounced  change  in  apparent  velocity  at  early  times.  These  results  suggest  that  all  GPR 
arrivals  are  sensitive  to  infiltration,  but  that  the  frequency  used  in  a  monitoring  experiment 
should  be  dependent  on  the  scale  of  heterogeneity  in  the  subsurface  and  the  length  of  time  of  the 
experiment  (i.e.,  maximum  depth  of  the  wetting  front). 

3.2.3  Numerical  Validation  of  Experimental  Results 

To  evaluate  the  empirical  data,  simulations  of  the  GPR  response  were  performed  using 
the  flow  model  HYDRUS-ID  [Simunek  et  ah,  2008]  and  a  MATLAB  based  2D  GPR  simulator 
created  by  Irving  and  Knight  [2006].  The  flow  was  simulated  using  the  known  boundary 
conditions  of  the  experiment,  i.e.,  the  applied  fluxes  given  in  Table  7  at  the  upper  boundary  and  a 
seepage  face  condition  on  the  lower  boundary.  The  independently  measured  hydraulic  properties 
of  the  sand  in  Table  6  were  used  in  the  simulations. 

A  comparison  of  the  observed  and  simulated  GPR  response  is  given  in  Figure  26.  The 
match  between  the  simulated  GPR  response  and  the  observed  data  is  excellent  given  that  the 
model  prediction  is  based  on  independent  samples  and  not  calibrated  to  the  observed  response. 
The  ground  wave,  second  wetting  front,  and  bottom  reflection  are  all  clearly  visible  in  the  data. 
An  exception  is  the  apparent  absence  of  the  reflection  from  the  first  wetting  front.  This  arrival 
can  be  found  in  the  simulated  data,  but  has  much  lower  amplitude  than  in  the  empirical  response 
so  is  difficult  to  identify  in  the  Figure.  The  patterns  observed  for  the  shifts  in  the  ground  wave 
arrival  also  show  excellent  qualitative  agreement  with  both  the  observed  and  simulated  average 
water  content  between  4-7cm  depth;  note  that  the  average  water  content  was  used  to  represent 
the  support  volume  of  the  capacitance  sensors.  The  traveltime  for  each  of  the  main  arrivals  in 
the  simulated  data  are  plotted  overtop  of  the  empirical  data  in  Figure  26a  to  provide  a  direct 
comparison  between  the  results.  It  is  apparent  that  there  is  generally  a  good  match  except  during 
the  period  between  50-75  minutes  where  the  simulated  traveltimes  are  less  than  the  observed 
response.  The  discrepancy  between  the  observed  and  simulated  water  content  can  be  explained 
by  the  underestimation  of  observed  water  content  by  the  simulations  in  the  near  surface  during 
this  period  (Figure  26c). 

An  important  outcome  of  the  empirical  and  numerical  results  is  the  observation  of 
hydrologic  trajectories  in  continuously  monitored  GPR  data  during  infiltration  events.  For 
example,  the  ground  wave  arrival  follows  a  complex,  but  well  defined  path  throughout  the 
experiment  that  can  be  reproduced  using  existing  numerical  models.  The  ability  to  reproduce 
complex  patterns  in  the  GPR  data  through  time  suggests  that  such  patterns  may  be  useful  tools 
for  hydrologic  characterization. 
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Observed  Data 


Simlulated  Data 


Mean  Water  Content  Between  4 -7 cm  Depth 


Figure  26:  (a)  Observed  GPR  response  during 
infiltration  event,  dotted  line  shows  comparison  to 
arrival  times  in  simulated  data;  (b)  simulated  GPR 
response  based  on  lab  measurements  of  hydraulic 
parameters,  (c)  average  observed  and  simulated 
water  content  between  4-7cm  depth,  (d)  water 
application  schedule. 


3.2.4  Coherency  analysis  of  GPR  data  using  hydrologic  trajectories 

The  observed  temporal  trajeetories  observed  in  the  GPR  response  to  infdtration  are 
reminiseent  of  the  spatial  trajeetories  observed  in  multi-offset  data.  Multi-offset  surveying,  sueh 
as  traditional  eentral  midpoint  (CMP)  and  wide  angle  reflection  and  refraction  (WARR)  surveys, 
are  well  known  methods  for  estimating  water  content  variations  with  depth  [Greaves  et  ah,  1996; 
van  Overmeeren  et  ah,  1997;  Garambois  et  ah,  2002],  Processing  of  multi-offset  data  is  typically 
performed  by  applying  normal  moveout  corrections  to  determine  the  one  dimensional  velocity 
structure  of  the  subsurface  [e.g.,  Yilmaz,  1987].  This  is  possible  because  the  change  in  reflector 
traveltime  with  antenna  offset  in  a  common  midpoint  gather  follows  a  hyperbolic  trajectory  that 
is  directly  dependent  on  RMS  velocity  (Figure  2).  Neidell  and  Taner  [1971]  discussed  the  use  of 
coherency  measures,  such  as  semblance,  calculated  along  these  moveout  trajectories  as  a  means 
for  determining  velocity  in  seismic  data  and  the  approach  has  since  also  become  a  standard 
method  for  velocity  estimation  in  GPR  workflows  [Fisher  et  ah,  1992;  Neal,  2004].  Water 
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content  can  be  calculated  from  interval  velocities  in  environments  with  low  electrical 
conductivity  when  an  appropriate  rock  physics  relationship  is  available  [Greaves  et  al.,  1996], 

Analogous  to  the  normal  moveout  approach  to  estimating  velocity,  I  suggest  that  the 
temporal  trajectories  in  constant  offset  GPR  data  can  be  used  to  identify  unsaturated  flow 
parameters.  Just  as  in  standard  semblance  analysis  to  estimate  velocity,  I  propose  that  it  is 
possible  to  calculate  the  coherency  of  the  GPR  signal  along  hydrologic  trajectories  that  are 
dependent  on  hydrologic  variables.  This  approach  to  analyzing  GPR  data  would  provide  a 
unique  means  for  identifying  these  variables  and  would  be  a  particularly  valuable  tool  if  it  could 
be  performed  by  defining  the  trajectories  using  methods  with  low  computational  cost. 

Predicting  the  GPR  traveltimes  for  different  arrivals,  e.g.,  ground  wave,  reflections, 
refractions,  etc.,  is  straightforward  for  models  with  sharp  wetting  fronts.  In  this  case,  the  same 
type  of  ray  calculations  performed  in  Section  3.1  can  be  used  to  determine  the  trajectories  as 
function  of  the  hydrologic  parameters.  For  each  set  of  hydrologic  parameters,  a  window  of  data 
along  the  trace  is  extracted  and  realigned  so  that  a  measure  of  coherency  between  the  traces,  such 
as  semblance,  can  be  calculated  (Figure  27). 


Time  [min] 


Figure  27:  (a)  Synthetic  GPR  data 
for  the  infiltration  example  using 
the  Philip  model.  The  dashed  line 
shows  the  semblance  analysis 
trajectory  for  the  wetting  front 
reflection  calculated  using  ray 
theory. 

(b)  GPR  data  extracted  in  a  1.5ns 
window  along  the  wetting  front 
reflection  trajectory  and  realigned. 


To  explore  the  use  of  semblance  calculations  for  different  types  of  arrivals  (i.e.,  ground 
wave,  wetting  front  reflection,  and  wetting  front  refraction),  a  sensitivity  analysis  was  performed 
for  each  of  the  hydrologic  variables  of  the  Philips  infiltration  and  rectangular  drainage  model 
described  in  Section  3.1.2.  The  results  of  the  sensitivity  analysis  are  shown  in  Figure  28  for  the 
parameter  set  given  in  Table  1  and  synthetic  GPR  data  in  Figure  27.  The  results  show  that  when 
all  three  arrivals  are  used  in  the  calculation  the  semblance  reaches  a  well-defined  maximum  at 
the  true  parameter  values,  except  for  the  N  parameter  which  has  a  poorly  defined  maximum 
consistent  with  the  results  of  the  MCMC  analyses  in  Section  3.1.  Taken  individually,  however, 
the  semblance  maximum  was  sometimes  poorly  defined  or  was  biased  toward  incorrect  values. 
This  was  particularly  true  for  the  ground  wave  results.  This  outcome  is  not  consistent  with  the 
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MCMC  results  discussed  earlier  where  the  ground  wave  was  able  to  provide  good  constraint  on 
the  hydrologic  model  parameters.  This  indicates  that  the  sensitivity  of  the  semblance  approach  is 
slightly  different  than  calibration  by  traveltimes. 
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Figure  28:  Sensitivity  analysis  for  semblance  as  a  function  of  parameters  of  the  Philip  infiltration 
and  rectangular  distribution  model.  The  bottom  three  rows  show  results  for  the  ground  wave, 
reflected  wave  from  the  wetting  front,  and  refracted  wave,  respectively.  The  top  row  shows 
results  by  summing  the  semblance  from  all  three  arrivals.  When  all  three  arrivals  are  combined, 
the  maximum  semblance  occurs  at  the  true  parameter  values. 


In  cases  where  more  complex  flow  occurs,  as  in  the  multi-flux  infiltration  test  performed 
in  the  lab  (Figure  26),  more  complicated  hydrologic  models  must  be  used.  In  this  case,  the  lack 
of  a  well-defined  wetting  front  in  the  model  may  prevent  the  calculation  of  arrival  times  using 
the  approach  described  above.  An  alternate  method  was  developed  in  this  research  that  is  based 
only  on  analysis  of  simulated  water  contents  derived  from  a  numerical  model.  First,  the  water 
contents  simulated  over  time  are  converted  to  dielectric  constants.  The  dielectric  constants  are 
used  to  calculate  profiles  of  reflection  coefficients  and  velocity  variations  with  depth.  The 
reflection  coefficients  are  indicative  of  areas  where  reflections  are  likely  to  occur,  i.e.,  strong 
signal  returns  would  be  expected  in  the  GPR  data.  The  velocities  are  then  used  to  remap  the 
reflection  coefficients  from  depth  profiles  to  profiles  in  terms  of  GPR  arrival  times.  The  result  is 
then  filtered  to  select  only  the  regions  with  the  strongest  reflection  coefficients,  which  gives  a 
map  of  trajectories  (as  in  Figure  27).  The  trajectories  are  then  used  to  extract  windows  out 
windows  from  the  observed  GPR  data  for  coherency  analysis.  Note  that  this  approach  is  based 
on  reflection  coefficients  and  therefore  only  provides  trajectories  related  to  reflections.  It  is 
therefore  most  suitable  for  cases  where  hydrologic  changes  dominate  the  GPR  response.  A 
sensitivity  analysis  for  the  multi-flux  infiltration  experiment  (Figure  26)  is  given  in  Figure  29.  In 
this  case  the  total  signal  power  in  the  window  was  used  as  the  coherency  measure  instead  of 
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semblance.  The  results  indicate  that  excellent  constraint  of  the  parameters  for  the  Mualem-van 
Genuchten  model  for  the  unsaturated  flow  parameters  can  be  achieved,  except  for  the  n 
parameter. 
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Figure  29:  Sensitivity  of  total  signal  power  along  trajectories  of  GPR  data  as  a  function  of 
hydrologic  parameters  for  the  multi-flux  infiltration  test.  Note  the  maximum  power  occurs  at 
the  true  parameter  value  (shown  with  a  dashed  line). 


The  sensitivity  analyses  presented  here  suggest  that  the  coherency  approach  to  estimating 
hydrologic  parameters  has  excellent  potential  for  the  hydraulic  analysis  of  GPR  data. 

3.3  Field-Scale  Study  of  GPR  Response  to  Infiltration 

Field-scale  experiments  have  also  been  conducted  to  perform  an  exploratory  investigation 
of  GPR  response  to  infiltration  in  heterogeneous  soils.  In  these  experiments,  we  used  a  sprinkler 
system  to  irrigate  an  approximately  10m  x  5m  region  of  a  silty-sand  soil.  Soil  moisture  probes 
were  installed  at  multiple  depths  to  monitor  the  movement  of  the  wetting  front.  The  rate  of  water 
application  over  the  site  was  monitored  through  time  using  a  tipping  bucket  rain  gauge.  GPR 
was  profiled  across  the  site  at  multiple  times  during  the  infiltration  experiment.  The  GPR  was 
profiled  between  fixed  stations,  but  due  to  an 
inability  to  integrate  in  real-time  positioning 
information  with  our  current  GPR  system  and  the 
discontinuation  of  odometer  wheel  triggering 
devices  by  the  manufacturer  of  our  radar,  we 
were  unable  to  precisely  position  the  data.  We 
found  that  careful  manual  positioning  at  fixed 
stations  was  a  poor  solution  to  the  positioning 
problem  because  for  each  GPR  profile  to  be 
collected  the  sprinklers  had  to  be  stopped.  As  a 
faster  alternative,  we  attempted  to  consistently 
move  the  antennas  at  a  constant  velocity  across 
the  site.  This  approach  introduced  the  potential 
for  positioning  errors.  Despite  the  possible 


Figure  30:  Two  10m  transects  of  GPR  data  collected  during 
a  field  infiltration  experiment.  The  upper  image  shows  the 
data  collected  prior  to  the  start  of  the  irrigation  sprinklers. 
The  bottom  image  shows  the  same  transect  profiled  after  16 
minutes  of  water  application  to  the  site.  Note  the  shift  in  the 
groundwave  and  other  reflections  is  a  result  of  increased 
water  content  (decreased  wave  velocity)  of  the  soils. 
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positioning  errors,  Figure  4e  illustrates  the  shift  in  the  arrival  time  of  the  groundwave  and 
different  refleetions  aeross  the  site  as  a  result  of  the  infiltration  event.  In  some  areas,  there  is  a 
change  in  the  character  of  the  GPR  image  that  is  likely  to  be  related  to  the  infdtration  (i.e.,  shifts 
in  reflections  to  later  times).  These  changes  are  not  uniformly  distributed  over  the  site  and  it  is 
currently  unclear  whether  this  is  an  effect  resulting  from  heterogeneous  water  infdtration  or 
errors  in  positioning  the  antennas.  This  experience  underscores  the  need  for  improved 
positioning  capabilities  for  this  type  of  field  research. 

In  addition  to  the  collection  of  GPR  data,  approximately  30  soil  samples  were  collected 
during  the  project  to  characterize  soil  variability.  Because  these  soils  drain  at  relatively  low 
pressure,  we  investigated  a  method  for  characterizing  the  soil  retention  curve  suitable  for  use 
with  large  numbers  of  samples.  Each  sample  was  packed  in  a  ~50cm  column  and  then  placed 
vertically  with  one  end  in  a  water  reservoir.  The  samples  were  allowed  to  equilibrate  such  that  at 
equilibrium  the  flux  in  the  column  is  zero  and  the  pressure  distribution  is  equal  to  the  elevation 
along  the  column.  The  columns  were  then  sampled  along  their  length  and  the  water  content  was 
measured  for  each  sample,  yielding  water  content  measurements  at  known  pressures.  While  the 
technique  showed  promise,  problems  with  packing  and  sampling  the  columns  put  into  question 
whether  the  results  are  reliable.  Additional  analysis  and  repeat  tests  are  required  to  fully  evaluate 
the  results. 


4.  Project  Conclusions 

This  short  term  research  project  has  demonstrated  the  potential  of  GPR  for  constraining 
the  parameters  of  unsaturated  flow  models  using  the  coupled  approach  to  hydrogeophysical 
estimation.  This  is  particularly  true  when  GPR  data  contain  insufficient  information  to 
independently  constrain  the  hydrologic  model.  In  this  case,  the  coupled  approach  allows  the 
hydrologic  model  to  regularize  the  inversion  and  produce  reliable  estimates  of  hydrologic 
parameters.  Stochastic  calibration  of  models  using  traveltime  data  was  shown  to  be  effective  for 
simple  flow  scenarios.  Coherency  analysis  along  hydrologic  trajectories,  a  new  methodology 
that  was  developed  as  a  result  of  this  research,  shows  promise  for  identifying  hydrologic  model 
parameters  for  more  complicated  flow  scenarios.  The  behavior  of  the  ground  wave  during 
infdtration  was  one  particular  point  of  interest  in  this  research.  It  was  shown  that  the  ground 
wave  response  is  sensitive  to  flow  processes  and  can  be  used  to  calibrate  hydrologic  models. 
This  research  also  suggests,  however,  that  reflection  data  provide  additional  information  on  the 
rate  of  flow  and  can  be  readily  included  in  model  calibration  algorithms.  Therefore,  it  is 
recommended  that  all  GPR  arrivals  be  used  to  calibrate  hydrologic  models  when  possible. 


36 


